library(osmdata)
library(opentripplanner)
library(sf)
library(tidyverse)
library(ggthemes)
library(ggspatial)
library(viridis)
library(dplyr)

Introduction

Tstops <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\T_stops\\MBTA_Rapid_Transit.shp")
affordable_2020 <- st_read("C:\\Users\\emmac\\Documents\\GitHub\\emmacolley-vis\\affordable_2020\\massbuilds-shp-20201004-a44e92.shp") %>%
  filter(MUNICIPAL == "Boston") %>%
  st_transform(crs = "+proj=longlat +datum=WGS84") #project to wgs84


MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"
opq(bbox = 'Boston MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/boston_streets.osm')
path_otp <- otp_dl_jar("OTP")
## The OTP will be saved to OTP/otp.jar
path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024)
## 2020-10-07 01:09:58 Basic checks completed, building graph, this may take a few minutes
## The graph will be saved to C:/Users/emmac/Documents/GitHub/emmacolley-vis/OTP
## 2020-10-07 01:10:12 Graph built
otp_setup(otp = path_otp, dir = path_data, memory = 1024)
## 2020-10-07 01:10:13 OTP is loading and may take a while to be useable
## Router http://localhost:8080/otp/routers/default exists
## 2020-10-07 01:11:14 OTP is ready to use Go to localhost:8080 in your browser to view the OTP
# Connect to opentripplanner
otpcon <- otp_connect()
## Router http://localhost:8080/otp/routers/default exists

1) Creating Isochrones

iso_5min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020,
                mode = "WALK", cutoffSec = 300) %>%
                st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")
## Warning in otp_isochrone(otpcon = otpcon, fromPlace =
## affordable_2020, mode = "WALK", : Failed to get isochrone
## with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 300},"id":"fid-24cdf481_175016df537_-7fac"}]}
iso_5min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = affordable_2020, 
                mode = "CAR", cutoffSec = 300) %>%
                st_transform(crs = MA_state_plane) %>%
  mutate(mode = "drive")

iso_all_modes <- rbind(iso_5min_drive, iso_5min_walk)
right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = affordable_2020) +
  geom_sf(data = Tstops, color = "lightslategrey") +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_manual(values = c("skyblue", "royalblue"), name = "Area that is reachable within", labels = c("5 min by car", "5 min by foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")
## Loading required namespace: raster

2) Plotting Walking and Driving distances

I’ve labeled each housing development with their neighborhood to begin to understand differences between neighborhoods and commuting distances.

iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(drive))) +
  geom_point(color = "lightslategrey", size = 3) +
  geom_text(aes(label = c("Allston", "South Boston", "Dorchester", "Dorchester", "5", "6", "7","8","9","10")), hjust=-.3) +
  scale_x_continuous(name = "Area within a ten-minute walking distance\nof housing development\n(square km)",
            breaks = breaks <- seq(10000, 130000, by = 20000),
            labels = breaks / 1000000) +
  scale_y_continuous(name = 
            "Area within a five-minute driving distance\nof housing development\n(square km)",
            breaks = breaks <- seq(0, 700000, by = 100000),
            labels = breaks / 1000000) +
  theme_minimal()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

3) Assigning a Transit Score to Each Isochrone

st_as_sf(iso_all_modes)
## Simple feature collection with 19 features and 4 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 228593.8 ymin: 893417.4 xmax: 237829.8 ymax: 901383.4
## CRS:            +proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs
## First 10 features:
##    id time              fromPlace  mode                       geometry
## 1   1  300 42.3165603,-71.0966731 drive MULTIPOLYGON (((232967.4 89...
## 2   1  300 42.3319784,-71.0502934 drive MULTIPOLYGON (((237064 8977...
## 3   1  300 42.3277973,-71.0563240 drive MULTIPOLYGON (((236345.8 89...
## 4   1  300 42.3049051,-71.0784474 drive MULTIPOLYGON (((234392.1 89...
## 5   1  300       42.3358,-71.1094 drive MULTIPOLYGON (((232188 8987...
## 6   1  300   42.360071,-71.149450 drive MULTIPOLYGON (((228881.4 90...
## 7   1  300 42.2974755,-71.1155956 drive MULTIPOLYGON (((231478 8935...
## 8   1  300 42.3007916,-71.1119202 drive MULTIPOLYGON (((231631.9 89...
## 9   1  300 42.3441666,-71.0636718 drive MULTIPOLYGON (((235949.3 89...
## 10  1  300 42.3387177,-71.0697751 drive MULTIPOLYGON (((235453 8984...
iso_transformed <- st_transform(iso_all_modes, crs = "+proj=longlat +datum=WGS84")
iso_transformed <- iso_transformed %>%
  mutate(transit_score = lengths(st_covers(geometry, Tstops)))
## although coordinates are longitude/latitude, st_covers assumes that they are planar
ggplot(iso_transformed) +
  annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
  geom_sf(aes(fill=transit_score), alpha=.5) +
  theme_map() 
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

4) Which Isochrones overlap with a T stop?

iso_transformed <- iso_transformed %>%
  mutate(num_Tstops = lengths(st_covers(iso_transformed, Tstops))) %>%
  mutate(has_Tstops = num_Tstops > 0)
## although coordinates are longitude/latitude, st_covers assumes that they are planar
n_Tstops_iso_transformed <- sum(iso_transformed$has_Tstops)

n_Tstops_iso_transformed
## [1] 6
ggplot(iso_transformed) +
  annotation_map_tile(zoomin = 1, type = "cartolight", progress = "none") +
  geom_sf(data = iso_transformed,
          aes(fill = has_Tstops)) +
  scale_fill_manual(values = c("gray85", "darkblue"),
          name = "Boston Neighborhoods\nby presence of a community center", 
          labels = c("Without a T stop",
                     "With a T stop")) +
  theme_map()
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

otp_stop()
## [1] "SUCCESS: The process \"java.exe\" with PID 21540 has been terminated."